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ABSTRACT 

Finite  element  methods  to  compute  approximate  solutions  to  flow  problems  involving 
the  flows  of  viscoelastic  fluids  are  discussed.  The  primary  goals  of  such  investigations  are  at 
least  three:  First,  to  evaluate  the  predictions  of  the  many  proposed  constitutive  theories  for 
viscoelastic  fluids.  Second,  to  model  measurement  flows  in  various  rheological  measurement 
devices  in  order  to  quantify  the  deviation  of  the  actual  flow  from  the  flow  which  must  be 
presumed  to  interpret  the  measurement.  Third,  it  is  hoped  that  these  methods  will  prove 
sufficiently  robust  to  allow  the  simulation  of  idealised  polymer  processes  with  the  aim 
of  aiding  in  the  design  of  such  processes  and  the  required  apparatus.  The  focus  of  the 
current  research  of  the  author  and  a  growing  number  of  others  is  on  two-dimensional, 
isothermal,  steady  flows  of  incompressible  fluids.  While  these  restrictions  will  be  seen  to 
be  non-essential  in  theory,  even  the  simplest  calculations  of  non-viscometric  flow  solutions 
will  be  seen  to  require  a  high  degree  of  computational  complexity  in  practice.  Nevertheless, 
the  current  finite  element  procedures  seem  to  show  promise  in  the  continuing  endeavor  to 
understand  this  challenging  class  of  problems.  . 
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FINITE  ELEMENT  METHODS  FOR  VISCOELASTIC  FLOW 

David  S.  Maikus* 


INTRODUCTION 

The  idea  to  apply  finite  element  methods  to  viscoelastic  fluid  flow  problems,  as 
opposed  to  finite  difference  methods,  seems  to  have  stemmed  from  the  desire  to  model  a 
variety  of  industrial  polymer  processes,  which  were  often  characterized  by  complex  flow 
geometries.  In  the  early  work  of  Tanner,  Nickel,  and  Bilger  [  1  ] ,  one  gets  the  clear  idea  that 
the  hope  was  to  develop  a  rather  general-purpose  code  which  could  be  used  by  the  industrial 
rheologist  to  tackle  geometrically  complicated  flow  situations,  particularly  those  flows  like 
the  extrusion  of  a  jet  of  liquid  material,  in  which  the  complex  geometry  itself  was  unknown  a 
priori.  Finite  elements  seem  to  be  ideal  tools  for  such  modelling  because  of  their  geometric 
flexibility.  Since  ref.  1  appeared,  it  has  been  discovered  that  there  were  a  number  of 
unanswered  questions  about  the  modelling  of  viscoelastic  fluids,  which  are  even  more  basic 
than  the  choice  of  discretization  scheme.  Many  of  the  problems  currently  being  solved  by 
this  author  and  others  are  of  an  idealized  nature  and  involve  little  geometric  complexity,  so 
that  they  could  just  as  easily  be  attacked  by  finite  differences,  and  the  questions  to  which 
we  seek  answers  could  as  well  be  asked  about  difference  methods  or  element  methods. 
There  are  a  number  of  investigators  who  have,  in  good  fluid-mechanical  tradition,  opted 
for  finite  difference  methods  [2-4].  Much  of  what  is  said  here  could  equally  apply  to  the 
difference  approach,  but  in  order  to  limit  the  present  discussion  to  the  manageable,  only 
finite  element  methods  will  be  discussed  here. 

The  reader  should  be  forewarned  that  the  author's  own  predilection  is  for  the  devel¬ 
opment  of  genera)  purpose  methods  which  can  be  applied  to  broad  classes  of  viscoelastic 
fluid  flow  problems  in  complex  geometries  without  reprogramming.  His  hope  is  that  once 
the  basic  difficulties  involving  choice  of  constitutive  equation,  convergence  of  nonlinear  it¬ 
eration  methods,  the  effect  of  singularities  of  unkown  order  —  and  others  still  —  are  more 
firmly  in  hand,  the  vision  on  ref.  1  may  be  more  completely  realized.  In  the  meantime, 
the  author  and  other  researchers  have  tended  to  solve  problems  which  do  not  tax  the  ge¬ 
ometric  flexibility  of  the  finite  element  formulation  (and  consequentially,  do  not  illustrate 
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its  power).  These  are  flows  which  might  be  classified  as  “measurement  flows”  in  many 
cases  —  very  well-characterised  flows  which  can  be  established  in  laboratory  situations 
for  the  purposes  of  making  measurements  of  material  properties.  The  lack  of  geometric 
complexity  leads  to  flows  which  can  be  idealized  as  flows  for  which  the  equations  of  motion 
can  be  solved  analytically;  however,  the  real  flow  is  some  perturbation  of  the  idealized  flow: 
Cone-and-plate  (5{  flow  actually  has  free  surfaces,  though  it  is  observed  that  this  does  not 
lead  to  an  unacceptable  disturbance  in  the  assumed  form  of  the  flow-field,  at  low  enough 
shear-rate.  The  entrance  into  a  slit-die  |5]  perturbs  the  flow-field  in  the  entry  region,  but 
the  velocity  field  and  pressure  gradient  soon  readjust  themselves  to  those  characteristic  of 
channel  flow,  at  least  at  low  enough  shear-rate.  There  are  many  other  such  flows,  includ¬ 
ing  plane  flow  over  a  transverse  slot,  discussed  below.  The  common  characteristic  of  such 
flows  is  that  they  are  perturbations  of  known  flows;  we  know  rather  well  what  qualitative 
features  the  numerical  results  must  display.  On  the  other  hand,  the  effect  the  perturbation 
of  the  assumed  flow  may  have  on  the  measurements,  particularly  at  the  extreme  range  of 
shear-rate,  is  a  question  which  seems  to  be  appropriately  addressed  by  numerical  mod¬ 
elling.  So  such  flows  seem  to  be  modest  testing  grounds  for  the  methods  under  discussion 
here,  yet  also  seem  to  raise  questions  to  which  numerical  modelling  may  be  able  to  supply 
useful  answers. 

We  will  see  that  investigations  of  such  problems,  humble  though  they  be  compared  to 
what  might  be  envisaged,  raise  thorny  questions  which  are  still  only  partially  resolved.  The 
questions  seem  to  have  at  their  heart  the  question  of  what  constitutive  equation  should  be 
used  and  why.  The  choice  of  constitutive  equation  will  be  seen  to  govern  a  basic  decision 
about  the  kind  of  finite  element  method  which  must  be  used.  It  will  be  observed  that 
different  constitutive  equations  can  lead  to  radically  different  behavior  of  the  numerical 
method,  yet  there  seems  to  be  little  to  differentiate  between  them  on  physical  grounds. 

CONSTITUTIVE  EQUATIONS 

Perhaps  the  most  basic  choice  in  the  design  of  a  numerical  method  devolves  from 
the  fact  that  there  are  two  seemingly  different  classes  of  constitutive  equations  which  have 
been  proposed  by  rheologists  and  which  seem  to  be  numerically  tractable:  the  differential 
and  the  integral  (6-9].  What  is  meant  specifically  by  these  terms  will  be  defined  shortly; 
the  important  thing  to  note  is  that  these  are  intersecting  classes  of  constitutive  equations, 
yet  neither  class  seems  to  entirely  contain  the  other.  At  first  sight,  the  two  classes  of 


constitutive  equations  seem  to  demand  totally  different  numerical  methods,  forcing  the 
numerical  modeller  to  make  an  irrevocable  choice  of  methodology  at  the  outset.  An  im¬ 
portant  point  to  be  observed  here  is  that  there  is  more  in  common  between  the  two  classes 
of  equations  than  first  meets  the  eye,  and  thus  the  potential  exists  for  a  more  unified 
numerical  approach  than  is  the  current  practice. 

General  Differential  Constitutive  Equation 

o  -  -pi  +  2(i{0)Re  +  (1  -  R)E 

_  v  (1) 

£  =  2^Tk 
k=y 

where  /i(0)  is  the  zero-shear  viscosity,  e  the  strain-rate  tensor,  R  <  1  is  a  retardation  ratio 
defined  by  R  =  A/T  for  a  characteristic  relaxation  time,  T,  and  retardation  time,  A  <  T, 
and  p  is  the  hydrostatic  pressure  function.  The  terms  on  the  right-hand  side  of  eq.  (1)  may 
be  interpreted  as  an  isotropic  contribution  to  the  stress,  the  Newtonian  contribution  to  the 
stress  —  possibly  arising  from  a  Newtonian  solvent  —  and  the  viscoelastic  contribution  to 
the  stress,  respectively. 

The  “extra  stress”  tensor,  E,  determines  the  non-Newtonian  contribution  to  the 
total  stress  tensor.  In  a  differential  constitutive  equation,  it  is  determined  by  one  auxiliary 
equation  for  each  of  the  “partial  stresses,”  r*  [8j, 

+  Bk  =  0 

At  dr  <2> 

_  =  +  (v .  V)t  -  r{Vv)T  -  (Vv)r 

where  A  *  is  the  kih  relaxation  time  in  a  relaxation  spectrum  (which  may  be  infinite). 
The  T  mentioned  above  can  be  taken  to  be  a  mean  relaxation  time  in  the  spectrum.  Bk 
is  a  tensor  function  which  may  depend  on  velocity  gradients  and  stresses(in  a  possibly 
nonlinear  fashion).  Specific  examples  will  be  cited  below.  Note  that  the  upper  convected 
time  derivative,  ^ ,  has  unknown  stresses  multiplying  unknown  velocity  derivatives,  and 
problems  with  such  constitutive  equations  are  inherently  nonlinear,  even  for  the  simplest, 
linear  choices  of  Bk.  With  the  obvious  correspondence  of  notations,  eqs.  (1)  and  (2)  are 
in  the  same  form  as  the  corresponding  equations  in  ref.  6,  generalized  slightly  in  order  to 
allow  the  possibility  of  more  than  one  partial  stress. 


Genera/  Single-Integral  Constitutive  Equation 

The  single-integral  designation  is  intended  to  make  a  distinction  between  the  consti¬ 
tutive  equations  under  consideration  here  and  those  involving  multiple,  iterated  integrals. 
In  fact,  the  theories  described  here  may  involve  sums  of  more  than  one  “single  integral.” 
While  it  makes  perfect  sense  to  consider  transient  flows  with  single-integral  constitutive 
equations,  the  computational  complexity  implied  by  such  flows  is  such  that  a  new  genera¬ 
tion  of  bigger,  faster  computers  would  probably  be  required.  Here  we  consider  steady  flows 
only,  whereas  no  such  restriction  was  imposed  on  our  discussion  of  differential  constitutive 
equations.  The  assumed  general  form  of  the  stress  field  is  given  by 

0  m  -pi  +  2 n{0)R4  +  (1  -  RW  (3) 

where  the  symbols  in  the  first  two  terms  are  the  same  as  those  in  eq.  (1).  The  extra  stress 
in  eq.  (3)  is  given  by  #*,  which  plays  an  analogous  role  to  £  in  eq.  (1).  For  the  integral 
constitutive  equations,  however,  the  extra  stress  is  given  by  an  explicit  calculating  formula, 
rather  than  implicitly,  as  in  eq.  (1).  It  is  assumed  to  have  the  form 

•'«£/  S  o(t)(r)mi(r)dT  (4) 

/=!  J~°° 

where  to/  is  a  memory  function,  and  So^  is  a  strain-measure.  The  memory  function  m;(r) 
has  the  following  form  in  most  examples  of  which  we  are  aware: 

m,(r)  =  T~x  ]T G^p^exvC^-)  (5) 

*= l 

The  constants  G fP  and  p^  whether  chosen  experimentally  or  theoretically,  must  be  chosen 
so  that  /^0O9(r)m/(r)dr  <  oo,  for  all  polynomials,  q,  of  arbitrary  degree.  This  is  automatic 
unless  Ni  =  oo,  which  does  occur  in  practice.  In  fact,  one  can  see  that  by  interchanging 
the  order  of  integration  in  eq.  (4)  with  the  summation  in  eq.  (5),  the  stress  contribution 
of  each  of  the  single  integrals  in  eq.  (4)  becomes  a  sum  of  partial  stresses,  each  associated 
with  a  member  of  the  relaxation  spectrum,  T/p*,  each  given  by  integration  of  a  strain- 
measure  against  a  single  exponential.  Thus  T/p*  in  eq.  (4)  and  A*  in  eq.  (2)  are  analogous 
relaxtion  times,  and  £  and  of  are  analogously  expressed  as  sums  of  partial  stresses. 

For  some  differential  constitutive  equations,  the  relation  between  them  and  a  re¬ 
arrangement  of  eq.  (4)  into  partial  stress  form  is  more  than  an  analogy:  the  integral 


constitutive  equation  amounts  to  a  direct  solution  of  the  hyperbolic  equations  (2)  for  a 
fixed  velocity  field  along  the  streamlines,  which  turn  out  to  be  a  double  family  of  real 
characteristics  [11].  Examples  of  this  will  be  pointed  out  subsequently.  It  is  important  to 
note  that  this  correspondence  between  integral  and  differentia)  constitutive  equations  only 
holds  in  very  special  cases.  There  are  constitutive  equations  derived  on  physical  grounds, 
which  have  an  integral  form,  yet  cannot  be  identified  as  characteristic  solutions  of  any  dif¬ 
ferential  constitutive  equation  given  by  eq.  (2).  Likewise,  there  are  differentia]  constitutive 

*  ;  ,  *  .  •  .  V  <  '  •-  * 

equations  which  do  not  seem  to  yield  to  closed-form  solutions  along  the  characteristics  in 
the  form  of  eq.  (4). 

The  strain  measures  in  eq.  (4)  are  assumed  to  have  the  following  dependencies: 

Sol')  =  S„<,>(E„(r),Eo'1M.^v|x(r)],Vy[x(0)])  (6) 

x(0)  is  the  position  at  which  o'  is  evaluated  in  the  current  configuration  of  the  fluid.  x(r) 
is  that  particle’s  position  at  time  r  in  the  past,  and  v  is  the  velocity  fieid  carrying  x(r) 
to  x(0),  inducing  the  stress  field.  Eo(r)  is  a  deformation-gradient  measuring  the  strain  in 
the  deformation  at  time  r  relative  to  time  0,  induced  by  v.  Eo(r)  may  be  the  usual  v 
or  something  similar  but  somewhat  different.  It  can  be  chosen  to  correspond  to  various 
convected  derivatives  in  differential  models  [8,9]  and  can  physically  represent  the  strain 

i 

due  to  molecules  which  undergo  motion  relative  to  the  motion  of  the  continuum  (non-afline 
motion)  [9,12]. 

To  compute  the  deformation-gradient,  it  is  neccessary  to  construct  the  path  x(r), 
which  is  a  streamline;  the  gradient  is  a  measure  of  strain  at  the  historical  time  and  position 
along  that  path  and  the  present  position.  The  following  equations  describe  the  path  and 
strain: 

x(r)  =  V|x(r)] 
x(0)  =  x0 

•  (7) 

E0(r)  =  F(x(r);  Vvjx(r)])E0(r) 

Eo(0)  =  I 

where  trF  =  0  for  all  r  <  0.  Note  that  this  is  a  system  of  equations  to  be  solved  in  reverse 
time.  The  first  group  is  possibly  nonlinear  and  determines  the  streamline;  the  second  group 
is  a  linear  system  whose  non-constant  coefficient  is  determined  by  the  streamline  solution. 
It  is  also  not  hard  to  deduce  that  Eo(r)  is  the  identity  in  the  present  configuration  and 
maintains  det(Eo(r))  =  1  for  all  historical  time  [10).  ,  < 


■* 
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Examples: 

1.  Differential  only: 

(a)  Phan-Tien /Tanner  (PTT)  [6] 

A*  1,  T  =  A,  ri  =  (1  -  R)p( 0) 

B  =■  Ba  -  rexp{  -  Irr }  -  2qe  +  A£(er  +  re) 

(  is  a  non-dimensional  parameter  of  the  model,  as  is  t;  When  the  latter  parameter 
is  not  zero,  there  is  no  known  integral  form  of  the  PTT  equation. 

(b)  Leonov  [6],  White-Metzner,  etc.  [7]. 

2.  Integral  only: 

(a)  Doi-Ed wards  [13] 

Af  =  1,  N  =  oo,  Pk  =  (2k  -  l)2 
Gk  =  C/Tpk,  F=Vv,  C=^l 
So  (r)  =  ^,(7, 7/)Col(r)  +  *2(/,//)C0(r) 

C0(r)  =  E0T(r)E0(r) 

The  scalar  functions  4>  j  and  <j>2  are  functions  of  the  two  nontrivial  principal  invariants, 
/  and  //,  of  Co*(r). 

(b)  Curtiss-Bird  [14] 

M  =  2,  A,  =  A,  =  oo,  =  p<2)  =  pk  =  (2k  -  l)2 

=  C\/Tpk,  G<2)  =  C2/p£,  F  =  Vv 
96p(0)  _  <192p(0) 

1  ir4(l  +  2c/3)  ’  2  ~  w2(l  +  2c/3) 

So(I^  is  the  same  as  So  of  (2a),  and  So(2^  is  a  more  complicated  tensor  function  of 
Co,  Cq1,  and  Vv(0)  1 10,14] .  <  is  the  “link-tension”  coefficient,  a  non-dimensional 
parameter  of  the  model;  with  c  =  0,  the  Curtiss-Bird  equation  reduces  to  the  Doi- 
Ed  wards. 

3.  Both  differential  and  integral: 

(a)  Oldroyd/ Jeffreys  (taking  A  =  1  for  simplicity,  not  necessity)  [8,9] 


T  =  A,  N  =  1 

B  =  r  -  2(1  -  R)p( 0)e  +  Aa(er  +  re) 


OTIC 


Note  that  taking  i-  a  and  <  =  0  in  (la)  yields  (3a). 

I  Copy  I  — . 

(b)  Johnson-Segalman  (with  a  nonzero  R)  (9,12]  X^specteo/1 

71  1  r  /.  ’  ;  '■  •  , 

T  =  A,  N  =  Af  =  1,  pi  =  1  H  . 

Gj  - /x(0)/A,  F  =  ae  -  u> 

So(t)  =  Ei1(r)EiT(r)  -  I 

(3b)  is  the  integral  form  of  (3a)  (and  thus  of  a  special  case  of  (la))  |9j.  u>  is  the 
vorticity  tensor.  The  parameter  a  is  a  “slip  parameter”  in  the  derivation  of  (3a), 
which  allows  for  non-afline  motion  of  a  strained  molecular  lattice  with  respect  to  the 
continuum  [9,12]. 

CHARACTERISTICS  AND  NUMERICAL  METHODS 

Steady  Plane  Flow  of  a  Johnson-Segalman  Fluid 

The  author  has  performed  computations  with  the  fluid  described  by  (3a)  and  (3b). 
Since  his  computations  employ  the  integral  form  of  the  constitutive  equation,  he  refers  to 
the  model  as  Johnson-Segalman.  It  turns  out  that,  for  the  purposes  of  analyzing  the  system 
of  equations  to  be  solved,  it  is  easier  to  deal  with  the  differential  form  (3b).  The  analysis 
discussed  here  has  only  been  carried  out  for  this  one  differentia]  constitutive  equation  (it 
does  include  several  simple  Maxwell  -  type  equations  for  particular  choices  of  a  [9]).  It  is 
hoped  that  the  conclusions  hold  for  other  differential  models  —  and  because  of  the  analogy 
alluded  to  earlier  —  to  integral  models  as  well.  In  what  follows,  take  R  =  0  and  N  —  1. 

The  assumed  form  of  the  equations  of  steady  motion  is 

V  •  a  +  f  =  p(u  ■  V)u  (8) 

where  u  is  the  (two-dimensional)  velocity  field,  a  is  the  stress  field,  f  a  body  force,  and 
p  the  fluid  density.  At  present  it  is  assumed  that  the  flow  is  a  plane  flow,  and  that  the 
material  is  exactly  incompressible, 

V  •  u  =  0  (9) 

For  the  Johnson-Segalman  fluid  with  parameters  as  specified,  the  equations  of  stress,  mo- 


tion,  and  continuity  take  the  form 


A^  +  B  =  0 
At 

p(u  •  V)u  +  Vp  -  V  ■  r  -  0 
V  u  =  0 


(10) 


The  analysis  of  ref.  11  begins  by  observing  that  in  this  differential  form,  eqs.  (10)  are  a 
quasilinear  system,  whose  characteristics  may  be  given  explicitly  in  terms  of  a  solution. 
The  analysis  of  ref.  13  reveals  that  below  a  critical  stress 

I.  For  a  fixed  u,  the  stress  equations  are  hyperbolic. 

II.  Streamlines  are  a  double  family  of  real  characteristics  of  the  stress  equations. 

III.  Solution  for  r  by  the  method  of  characteristics  yields  a  problem  which  is  elliptic  in  u. 
The  stress  is  given  by  an  integral  stress  calculator  of  the  form  of  eq.  (4). 

It  should  be  noted  that  the  sense  in  which  the  equations  of  motion  are  elliptic  must  be 
carefully  qualified:  There  are  no  real  characteristics  of  the  equations  of  motion  alone,  be¬ 
low  the  “sonic  transition”  (critical  stress).  But  since  there  may  be  characteristics  entering 
through  boundaries  at  inflows  (with  corresponding  outflows),  stress  data  —  or  the  equiva¬ 
lent  —  must  be  somehow  specified  there.  So  the  whole  problem  is  not  elliptic  in  the  usual 
sense,  when  characteristics  cross  domain  boundaries. 

There  is  a  way  to  make  the  equations  of  motion  elliptic  in  the  usual  sense  when  a 
given  stress  field  is  specified,  ajtfi  that  is  to  rewrite  eq.  (10)  in  equivalent  form, 

V 

a£T  +  b  =  o 

At 

p(u- V)u+ Vp-pV2u  =  Vt-/*V2u  (10a) 

V  •  u  =  0 


in  which  p  is  a  positive,  viscosity-like  constant.  In  an  iterative  solution  method  which 
uses  the  right-hand  side  of  the  middle  equations  of  eqs.  (10a)  from  a  previous  velocity 
iterate  and  the  stresses  calculated  from  the  hyperbolic  stress  equations  in  that  velocity 
field,  solution  for  the  next  velocity  iterate  can  be  a  truly  elliptic  problem.  This  is  exploited 
in  some  of  the  methods  to  which  we  now  turn  our  attention. 
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The  Spectrum  of  Discrete  Methods 


The  results  of  ref.  1 1  have  interesting  ramifications  for  numerical  methods  for  vis¬ 
coelastic  fluids.  In  what  follows,  we  assume  that  computations  are  to  be  carried  out  in  a 
steady  flow  regime  below  the  critical  stress  for  type-change  (numerical  experiments  by  the 
author  suggest  that  with  the  contitutive  equations  currently  employed,  the  critical  stress 
is  so  large  as  to  be  essentially  inaccessible  by  current  numerical  methods).  We  assume  that 
the  characteristics  of  the  equations  to  be  solved  for  any  constitutive  equation  of  the  forms 
considered  in  this  paper  are  the  same  as  they  are  in  the  Johnson-Segalman  case.  While 
the  analysis  has  not  been  carried  out  in  such  generality,  this  assumption  does  not  seem 
unwarranted  in  the  light  of  numerical  evidence.  The  methods  currently  employed  by  most 
researchers  fall  somewhere  on  the  spectrum  of  methods  given  in  Table  1. 


METHOD 

MIXED 

STRESS  ALT.  v,p 

PURE  v,p 

KEY 

FEATURE 

AU  *1  once 

Variant  I  J  Variant  2 

Total  numerical  j  Metli.  chars. 

PDE 

TYPE 

Mixed 

1 

Hyperbolic/Elliptic  J  Hyperbolic/Elliptic 

UNK8. 

u.r.p 

1 

1 

r  1  r 

u>P  ]  u,p 

1 

1 

1 

i 

U>P  j  u 

1 

1 

REFS. 

m 

i 

1 

[is)  ;  [is] 

i 

1 

_ 1 - 

m 

Table  1 

Spectrum  of  Finite  Element  Methods  for  Viscoelastic  Flow. 


t  the  extreme  left  end.  there  is  a  class  of  methods  applicable,  in  principle,  to  any  consti- 
itive  equation  which  has  a  differential  form.  These  methods  treat  the  quasilinear  system 
i  a  system  of  mixed  type  and  discretize  velocity,  pressure,  and  extra  stress  independently 
a  mixed/Galerkin  method.  At  the  extreme  right  end  of  the  spectrum  are  methods  which 
>ply  to  those  constitutive  equations  which  have  an  integral  form.  They  discretize  only 
ial  velocity  fields  in  an  iterative  sequence  of  approximations  to  the  solution  of  the  equa- 
ons  of  motion  and  continuity;  even  the  pressure  can  be  eliminated  if  a  penalty  approach  is 
ted.  In  these  right-wing  methods,  stress  is  computed  via  the  integral  stress  calculator  of 
|.  (4)  as  a  functional  of  the  trial  velocity  fields.  Each  of  the  extreme  right-  and  left-wing 
ethods  have  their  own  distinct  advantages  and  disadvantages.  Left-wing,  mixed  methods 
*d  many  stress  unknowns  —  particularly  if  more  than  one  relaxation  time  and  corre- 
tonding  partial  stress  must  be  discretized.  Left-wing  methods  thus  involve  large  matrices 
ith  large  bandwidths.  but  transient  problems  can  be  solved  by  well-known  time-stepping 
ethods.  Left-wing  methods  also  require  the  specification  of  boundary  data  where  charac- 
ristics  emanate  from  boundaries.  How  to  specify  such  data  for  partial  stresses  is  a  subtle 
atter  [11],  and  the  possiblity  of  over  or  under  specifying  boundary  data  is  real  one.  On 
ie  other  hand,  left-wing  methods  give  the  Galerkin  residual  in  terms  of  stress  unknowns, 
hich  are  primary  variables;  thus  differentiation  of  the  Galerkin  residuals  with  respect  to 
>dal  variables  as  required  to  form  Newton’s  method  for  solving  the  discrete  nonlinear 
alerkin  equations  is  straight-forward. 

The  right-wing  methods  make  some  direct  trades  of  advantage  for  disadvantage  when 
>m pared  to  the  left-wing  methods.  Only  velocities  need  be  discretized;  computational 
ark  is  not  increased  by  adding  more  relaxation  times.  The  specification  of  stress  on 
coming  characteristics  can  be  handled  implicitly,  in  a  way  that  seems  to  work  rather 
ill.  if  not  theoretically  justified  at  present.  The  incoming  streamlines  are  assumed  to  be 
ose  of  a  known  flow  extending  to  infinity  in  a  known  geometry  (such  as  channel  flow),  and 
us  the  required  contributions  to  the  integrand  of  eqs.  (4)  and  (11)  are  actually  known 
lalytically  once  V  *'-'njc|e  enters  this  “predecessor  flow”  j  10,19! .  On  the  other  hand, 
ansient  methods  seen  rly  impossible,  and  Newton's  method  seems  difficult,  at  best, 
cause  of  the  complicate  inctiona)  relation  between  the  stresses  in  the  Galerkin  residual 
d  the  velocity  field.  1  .nost  cases,  it  seems  that  the  choice  of  left-wing  and  right-wing 
ethod  is  made  a-prio  by  the  form  of  the  constitutive  equation:  integral  models  seem  to 
mand  right-wing  me  .hods  and  differential  left-wing  methods. 
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Figure  7 

Flow  with  the  Concocted  Constitutive  Equation  at  Df  =  3.7. 

The  flow  pictured  in  Figure  6  has  De  =  2.5,  L0  =  1.8,  and  viscosity  reduced  to  46%  of  its 
zero-shear  value.  The  flow  of  Figure  7  is  at  De  =  3.7.  but  Lo  =  1.6,  which  is  lower  than 
the  stress-ratio  of  Figure  6:  in  Figure  7,  the  viscosity  is  40%  of  its  zero-shear  value.  The 
stress-ratio  of  the  concocted  fluid  model  goes  through  a  maximum,  due  to  the  retardation 
Lime;  the  flow  of  Figure  6  has  shear-flow  at  the  wall  which  is  just  about  at  that  maximium. 
16  Laguerre  integration  points  were  used  for  the  quadrature  of  eq.  (ll)  when  computing 
;he  flows  of  Figures  6  and  7. 

In  both  the  flows  of  Figure  6  and  Figure  7.  the  Broyden  iteration  scheme  had  a 
ougher  time  converging  than  it  did  with  the  Curtiss-Bird  fluid  at  a  much  higher  De  and 
L0.  The  De  =  3.7  case  took  40  iterations  to  get  only  a  factor  of  25  reduction  in  the  residual. 
\ttempts  to  converge  at  higher  De  failed.  The  De  —  2.5  case  took  20  iterations  to  get  a 
actor  of  20  reduction  in  the  residual.  Thereare  some  qualitative  differences  between  the 
lows  of  Figures  6  and  7  and  the  flow  of  Figure  5.  These  may  be  related  to  the  greater 
lifficulty  in  convergence,  though  this  cannot  be  stated  with  any  certainty.  We  note  that 
he  streamlines  just  outside  the  slot  are.  if  anything,  more  distorted  by  the  presence  of  the 


damping  function 


<t>(II)  =  txp{-cvll  -  3)  (22) 

which  just  multiplies  the  strain-measure,  So,  of  the  Johnson-Segalman  model.  II  is  the 
second  principal  invariant  of  the  non-affine  strain  measure,  Eq  1E0  T  No  physical  reasoning 
lies  behind  this;  damping  functions  are  intended  to  account  for  the  effect  of  breaking 
temporary  junctions  in  the  lattice  under  large  strain.  To  what  extent  this  would  happen 
if  non-affine  motion  of  the  lattice  is  allowed,  or,  if  so,  whether  the  damping  function  for  a 
non-affine  strain-measure  should  take  the  form  of  eq.  (24)  is  unknown.  At  any  rate,  the 
flow  of  Figures  6  and  7  use  this  concocted  model,  and  take  c  =  0.1;  the  other  parameters 
were  as  just  described  for  the  Johnson-Segalman  with  a  retardation  time. 


constitutive  equations  [2,7,16,20j.  The  run  was  terminated  when  the  residua]  measure  of 
ref.  10  had  been  reduced  by  a  factor  of  84  (a  similar  reduction  factor  to  those  reported 
in  ref.  10).  The  residual  appeared  to  be  continuing  to  converge,  but  the  computation  was 
terminated  because  the  experience  of  ref.  10  shows  that  most  quantities  of  interest  have 
stabilized  to  a  couple  of  decimal  digits  with  a  comparable  residual  reduction.  This  residual 
reduction  was  accomplished  in  18  inverse  Broyden  iterations.  The  computation  took  15 
hours  of  VAX  11/780  CPU  time  on  a  VAX  with  a  standard  configuration,  using  double 
precision  hardware. 

There  are  several  interesting  features  of  Figure  5.  The  flow  is  from  right  to  left,  and 
we  observe  substantial  vortex  distortion  of  the  type  observed  in  ref.  4.  We  also  note  that 
there  is  little  disturbance  in  the  streamlines  near  the  mouth  of  the  slot  when  compared 
to  those  in  the  Stokes-flow  of  Figure  4  or  the  streamlines  in  the  flow  visualizations  of  ref. 
4.  Also,  of  the  five  equally-spaced  streamlines  in  the  slot,  three  are  crowded  together 
near  the  dividing  streamline;  this  indicates  that  the  velocity  of  fluid  in  the  vortex  is  lower 
compared  to  the  velocity  at  the  dividing  streamline  in  the  non-Newtonian  case  than  it  is 
in  Stokes-flow.  This  appears  to  be  due  to  shear-thinning;  the  viscosity  in  undisturbed  flow 
at  the  wall  in  the  case  of  Figure  5  is  only  9%  of  its  zero-shear  value.  The  very  high  driving 
stresses  in  the  die  are  not  transmitted  to  the  slot,  and  the  fluid  slides  by  without  feeling 
much  effect  of  the  slot. 

We  turn  our  attention  to  the  flow  of  another  fluid  in  the  same  geometry,  using  the 
same  mesh  as  before.  The  author  really  would  have  liked  to  have  shown  results  involving 
a  Johnson-Segalman  fluid  at  significant  De,  but  these  are  not  available.  Johnson  and 
Segalman  found  that  a  value  of  their  slip-parameter  of  a  =  0.87  gave  good  agreement  of 
viscometric  functions  with  published  data  [12].  With  this  value  of  a,  a  very  high  zero  shear 
viscosity  (implying  negligible  Reynolds  number),  and  R  =  0  in  eq.  (3),  the  inverse  Broyden 
iterations  of  eq.  (18)  diverge  dramatically  at  about  Dt  =  1.0.  Adding  a  substantial 
Newtonian  viscosity  with  R  —  1/3  helps  some,  as  one  might  expect,  but  divergence  ensues 
at  a  disappointing  De  =  1.5.  In  order  to  obtain  an  integral  constitutive  equation  which 
shear  thins  less  than  Curtiss-Bird,  behaves  better  than  Maxwell,  does  not  imply  very 
peculiar  lattice  motion  which  would  be  the  result  of  taking  a  very  small  slip-parameter  (9j, 
and  which  allows  computation  with  at  least  a  moderate  Dt.  the  author  just  tinkered  with 
the  Johnson-Segalman  equation.  Favorable  results  were  obtained  by  adding  a  Wagner-type 


De  is  the  “Deborah  number"  and  L0  the  “stress-ratio."  -y  is  the  shear-rate  at  the  wall  of  the 
die,  far  from  the  slot,  and  JVj  and  a  are  taken  at  the  same  location.  For  the  flow  pictured 
in  Figure  5,  Dt  —  17.9,  and  L0  =  4.1;  the  computation  used  the  Curtiss-Bird  equation 
1 10,14],  with  the  same  parameters  given  in  ref.  10  for  the  “EFN  Fluid."  This  is  intended 
to  model  a  high-viscosity,  concentrated  polymer  system  without  significant  branching  of 
the  molecular  chains.  The  resulting  viscosity  is  so  high  («  104  poise  at  zero  shear)  that 
the  Reynolds  number  is  negligible,  even  when  shear-thinning  is  accounted  for. 


The  computations  of  ref.  10  were  carried  out  with  Gaussian  memory  formulas  (eq.  (11)) 
with  ten  points;  computations  were  reported  there  up  to  Dt  =  7.5.  It  has  since  been 
determined  that  the  ten-point  Gaussian  formula  loses  its  accuracy  in  simple  shearing  flow 
(and  thus  in  the  undisturbed  die-flow)  for  Dt  much  larger  than  that.  Furthermore,  one 
observes  that  the  ten-point  formula  is  inaccurate  enough  to  lead  to  an  artificial  shear- 
stress  maximum  f 30,20]  for  very  high  Dt.  For  that  reason,  the  computations  leading  to 
Figure  5  were  carried  out  with  18-point  formulas  for  each  memory-integral.  With  the 
18-point  formula,  no  difficulty  in  convergence  of  the  Broyden  iterations  was  observed  — 
no  “barrier*  or  “high  Weissenberg  number  problem"  which  has  been  observed  for  other 


Stokes-flow  in  the  Domain  of  Figure  2,  Using  the  Mesh  of  Figure  3. 


We  note  that  the  symmmetry  of  the  streamlines  is  perfect  to  numerical  accuracy,  and  any 
very  slight  deviations  from  symmetry  in  Figure  4  arise  from  the  streamline  integration  used 
for  graphical  purposes  [25].  This  integration  does  not  make  use  of  the  element-wise  conic 
sections  on  triangles,  but  an  approximate  method.  In  all  streamline  plots  shown  here,  the 
streamlines  are  not  equally  spaced  contours  throughout  because  the  vortex  in  the  slot  is 
so  faint  in  comparison  to  the  flow  outside.  Instead,  there  are  ten  equally  spaced  contours 
outside  the  slot,  the  streamline  at  the  wall  in  which  the  slot  is  cut,  and  five  equally  spaced 
contours  in  the  slot.  These  latter  contours  have  values  which  evenly  divide  the  range  of 
values  between  the  streamline  at  the  wall  and  the  maximum  value  of  the  stream  function. 

We  turn  our  attention  to  the  non-Newtonian  flow  illustrated  in  Figure  5.  Among 
other  things,  Figure  5  illustrates  the  fact  that  there  are  some  constitutive  equations  with 
which  rather  high  shear-rates  can  be  obtained.  Two  non-dimensional  numbers  are  useful 
in  quantifying  the  degree  of  nonlinearity  due  to  non-Newtonian  effects, 

De=7(T-A) 
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Mesh  Discretizing  the  Domain  of  Figure  2  with  1008  Macroelements. 

The  distribution  of  elements  was  determined  from  extensive  numerical  experimentation, 
and  the  author  believes  it  produces  more  accurate  results  than  do  any  of  those  of  ref.  10. 
The  domain  parameters  of  Figure  2  are  h  =  1,  b  =  0.5,  and  q  =  2.  Figure  4  shows  a  plot 
of  the  streamlines  in  Stokes-flow  (Newtonian,  zero  Reynolds  number)  obtained  using  this 
mesh;  this  will  serve  as  a  useful  baseline  in  what  follows. 


't  A*  «- 


a  long  history  and  is  discussed  in  some  detail  in  refs.  3.4.10,17,20.  and  22*25.  It  is  of 
practical  importance,  because  it  seems  feasible  to  design  measurement  devices  based  on 
the  relation  between  Pt  and  Ni .  which  predict  the  latter  based  on  continuous  measurement 
of  the  former  and  the  wall  shear-stress  in  undisturbed  flow,  a  10,23,24).  It  has  been  found 
experimentally  that  the  P,  A’j  relation  accurately  describes  the  behavior  of  various 
viscoelastic  fluids  in  flows  over  slots  23,24..  This  appears  to  happen  in  spite  of  the  fact 
that  the  assumptions  behind  the  derivation  22  of  eq.  (20)  are  violated  ]  10,23].  The 
purpose  of  the  investigation  undertaken  by  the  author  is  to  explore  this  apparent  anomoiy 
in  detail,  using  his  numerical  techniques.  Significant  progress  has  been  made  in  these 
investigations  since  the  appearence  of  ref.  10;  however,  these  results  are  best  presented 
in  a  future  publication.  Here  we  will  focus  on  some  interesting  qualitative  features  of 
plane  flow  over  a  transverse  slot  and  the  apparent  relation  of  these  qualitative  features 
to  the  numerically  observed  hole-pressure.  The  results  will  also  illustrate  the  dramatic 
and  somewhat  puzzling  differences  resulting  from  using  different  constitutive  equations  to 
study  the  same  flow. 

Numer/ca/  Results 

Figure  3  illustrates  a  mesh  of  1008  crossed-triangle  macroelements,  which  has  been 
designed  by  the  author  to  discretize  the  domain  of  Figure  2  and  to  investigate  hole-pressure 
problem  of  the  previous  subsection. 
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from  which  one  may  deduce  that  if  Ha  were  the  inverse  Jacobian,  eqs.  (18)  would  be 
Newton’s  method.  The  inverse  Broyden  method  given  by  eqs.  (18)  is  a  fairly  robust 
method  for  computations  with  some  constitutive  equations,  but  it  diverges  at  rather  low 
Deborah  numbers  with  others.  The  reason  for  this  is  not  entirely  clear  at  this  point,  but 
one  avenue  worth  exploring  fully  is  the  possibility  of  finding  more  accurate  approximations 
to  the  actual  inverse  Jacobian,  and  the  author  is  currently  involved  in  studying  this  matter. 

A  MODEL  PROBLEM 

Plane  Flow  over  a  Transverse  Slot 


Figure  2 

Domain  for  Plane  Flow  over  a  Transverse  Slot. 

Here  we  consider  the  flow  of  a  viscoelastic  fluid  over  a  transverse  slot  {3,4,10,17,20,22- 
25].  The  author  [10,26]  has  been  investigating  the  prediction  that  at  negligible  Reynolds 
number  (based  on  slot  width),  the  difference  between  the  thrust  on  the  wall  at  the  top  and 
bottom  of  the  slot,  pictured  in  Figure  2 

P<  =  Pi-  Pi  (20) 

is  related  to  the  first  normal  stress  difference,  N\,  in  a  simple  fashion.  This  problem  has 
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Solving  the  Nonlinear  Equations 

Here  we  briefly  discuss  the  kind  of  numerical  method  which  can  be  used  to  find  the 
zero(s)  of  eq.  (14).  One  may  deduce  that  the  functional  dependence  of  the  right-hand 
side  of  eq.  (14)  on  the  unknown  velocity  field  is  so  complicated  that  Newton’s  method  is 
difficult  at  best.  On  the  other  hand,  it  has  been  found  that  iterative  methods  based  on 
the  finite  element  approximation  to  the  Stokes  operator,  which  puts  all  nonlinearity  on 
the  right-hand  side,  do  not  have  very  favorable  convergence  properties,  even  at  moderate 
Deborah  numbers  [2,7,19] .  The  method  illustrated  here  is  a  Gauss-Newton  type  which  be¬ 
gins  from  an  estimate  of  the  solution,  u0,  which  represents  nodal  values  of  a  finite  element 
solution  computed  using  the  Stokes  operator,  but  with  inflow  and  outflow  boundary  condi¬ 
tions  appropriate  to  the  non-Newtonian  solution  (the  “pseudo-Newtonian”  solution).  The 
iteration  matrix  is  an  approximation  to  the  unavailable  inverse  Jacobian.  The  inverse  Jaco¬ 
bian  approximation  begins  using  the  Stokes  approximation  and  continues,  using  rank-one 
updates,  ri  and  vj,  to  better  approximate  the  inverse  Jacobian  as  the  iteration  progresses. 
This  method  is  known  as  the  “inverse  Broyden”  method,  and  it  is  quite  efficient,  as  it  only 
requires  a  single  matrix  elimination  proceedure  to  be  carried  out  on  H0. 

for  k  =  0 

H0  =  inverse  Stokes’  matrix, 

Uq  =  pseudo-Newtonian  solution,  (17) 

R(u0)  =  residual  of  eq.  (15)  at  Uo 
s0  =  -H0R(u0) 

for  k  =  0, 1, 2, . . . 

Uk+i  =  uk  4  sk, 

yk  =  R(«k+i)  -  R(uk), 
vk  =  HkT8k/(skrHkyk), 
rk  =  sk  -  Hkyk, 

Hk.i-Ho  +  E*=oriV|7\ 

8k- 1  =  ll  ~  vkrR(uk4j)]rk, 
next  k. 

By  rearrangement  of  eq.  (18),  it  follows  that 

Uk-l  -  uk  =  8k 

•k  =  HkR(uk).  (19) 


The  Stream  and  Drift  Functions 


The  construction  of  the  particle  paths  of  Figure  1  and  the  accumulation  of  the 
deformation-gradient  via  eq.  (13)  are  the  dominant  computational  cost  of  the  method 
presented  here.  There  are  several  interesting  aspects  of  the  chosen  finite  elements  which 
contribute  to  the  efficiency  with  which  this  can  be  done.  Since  the  interpolating  functions 
determining  the  test  and  trial  functions  on  each  triangle  are  linear,  it  follows  that  there  is 
a  quadratic  stream  function  on  each  triangle,  t/>(x)  (it  actually  turns  out  to  be  a  quadratic 
spline  on  the  whole  mesh).  Thus  the  boundary  crossings  illustrated  in  Figure  1  can  be 
found  in  sequential  triangles  upstream  of  each  ( t  using  the  quadratic  formula, 

V»(xi)  =  0(x3) 

„  .  (15) 

l{xl2,x\)  -  0 

The  relation  involving  /(•,•)  is  the  requirement  that  X3  be  on  the  line  defining  one  of  the 
boundaries  of  the  current  triangle.  Thus  eq.  (15)  is  a  quadratic  equation  in  one  of  the 
components  of  X3,  either  xj  or  x\-  To  associate  a  historical  time  with  Xa,  it  turns  out  that 
for  the  chosen  elements,  the  time  of  transit  between  two  points  on  the  same  streamline 
can  be  determined  from  the  drift  function,  which  is  known  analytically  [10] . 

r,  -  r2  =  w(x j)  -  u>(x3)  (16) 

The  drift  function  is  denoted  by  w,  its  analytic  expression  involves  various  cases  depend¬ 
ing  on  the  characteristics  of  the  pathline  equations  in  eq.  (7)  [10].  Stress  calculation  can 
thus  proceed  by  finding  boundary  crossing  points  only  and  accumulating  strain  at  these 
points  via  eq.  (13),  using  the  times  associated  with  the  boundary  crossings  computed  via 
eq.  (16).  The  accumulation  proceeds  until  the  boundary  crossing  times  in  one  element, 
T\  and  T2,  bracket  one  or  more  temporal  integration  points  of  eq.  (11),  at  which  time(s) 
the  integrand  of  eq.  (11)  can  be  evaluated  using  the  analytic  expression  for  the  strain  1 10 
and  eq.  (13)  with  r  taken  as  the  value  of  the  temporal  integration  point(s).  Since  there 
are  potentially  six  solutions  to  eqs.  (15),  when  each  boundary  of  the  current  triangle  is 
considered  as  a  possible  candidate  on  which  to  find  x3.  the  drift  function  is  also  indispen- 
sible  in  determining  the  actual  boundary  crossing  of  a  particle  as  the  one  which  gives  the 
maximum  value  of  rj  satisfying  r2  <  ri . 


at  discrete  points  in  historical  time.  r*.  For  a  single  relaxation  time,  or  N  =  1  in  eq. 
(5),  the  Gaussian  formula  is  given  by  the  well-known  Laguerre  family  of  formulas.  For 
other  memory  functions  (assumed  to  be  sums  of  exponentials  here),  analogous  formulas 
can  be  derived  {10].  It  should  be  pointed  out  that  the  number,  Nt„  of  integration  points 
is  not  explicitly  related  to  the  number,  N.  of  relaxation  times  (which  is  often  infinite). 
This  is  in  distinct  contrast  to  the  situation  with  differential  constitutive  equations,  where 
each  relaxation  time  adds  a  new  partial  stress  and  a  new  set  of  stress  equations  to  be 
solved.  With  the  integral  form,  adding  relaxation  times  need  not  increase  Np  and  the 
corresponding  computational  effort.  Experience  shows,  however,  that  if  adding  relaxation 
times  significantly  spreads  out  the  relaxation  spectrum,  an  increased  Np  may  be  required 
to  maintain  accuracy. 


Evaluating  the  Strain 

The  author  has  chosen  a  basic  approach  to  solving  the  tracking  equations  (7)  which 
involves  the  choice  of  finite  elements  for  which  the  solution  for  the  particle  pathline  (stream¬ 
line)  is  known  exactly  in  piecewise  fashion  as  it  crosses  each  element.  For  this  purpose, 
the  quadrilateral  element  composed  of  four  linear  triangles  defining  the  diagonals  of  the 
quadrilateral  have  been  chosen;  these  are  illustrated  in  Figure  1.  The  composite  pattern  is 
chosen  to  produce  accurate  satisfaction  of  the  continuity  equation  with  an  element  of  suffi¬ 
cient  over-all  accuracy.  Linear  triangles  in  other  arrangements  will  not  do  this.  In  order  to 
compute  the  integrand  of  eq.  (11),  we  observe  that  if  a  pathline  is  known,  the  strain  mea¬ 
sures  of  eq.  (6)  are  readily  computed  from  nodal  variables  and  assumed  interpolations  if 
Eo(r)  is  known.  Eo  is  the  solution  to  the  evolution  equations  for  the  deformation-gradient, 
which  are  part  of  eqs.  (7).  This  solution  is  known  piecewise  analytically  on  each  element 
of  the  type  which  is  employed.  We  proceed  as  follows:  We  define  Eo(rj)  to  be  the  solution 
to  the  strain  evolution  equations  back  to  historical  time  t j,  at  which  the  particle  being 
tracked  is  located  at  the  boundary  of  some  triangle  in  the  mesh.  We  define  ETl(r)  to  be 
the  solution  to  the  following  strain  evolution  equation,  which  is  derived  from  eq.  (7)  by 
modifying  the  initial  condition  to  an  interface  condition: 


Er,  =  FEfl 
Er.(r,)  =  I 


(12) 


It  turns  out  that  Er,  can  be  easily  determined  analytically  for  the  chosen  elements  [10]. 
The  whole  deformation -gradient  back  to  r  <  rj  can  then  be  mutiplicatively  accumulated 
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fundamental  difference  between  integral  and  differential  forms  is  only  the  question  as  to 
whether  a  closed-form  quadrature  is  known  for  the  stresses  or  integration  of  ODEs  along 
the  streamlines  is  required. 

An  alternative  way  to  look  at  Table  1  is  as  an  indicator  of  the  degree  of  “La- 
grangianess”  of  the  various  methods.  From  that  point  of  view,  the  left  and  right  wings 
of  the  table  do  not  define  extremes  any  longer.  One  may  properly  consider  the  left-wing 
methods  to  be  totally  Eulerian  in  that  all  quantities  solved  for  are  referred  to  fixed  spa¬ 
tial  locations  and  not  to  particles,  as  we  move  rightward  on  the  table,  streamlines  — 
and  thus  particle  paths  —  begin  to  play  an  increasingly  important  role.  But  the  right- 
wing  methods  still  retain  the  equations  of  motion  in  Eulerian  form  and  are  only  partially 
Lagrangian  in  their  need  to  track  particles  in  order  to  compute  stresses.  At  the  far  right 
wing  of  Lagrangianess  is  found  Hassager’s  method  j21),  which  uses  a  deforming  mesh.  This 
method  has  some  drawbacks  —  particularly  in  flows  with  recirculations  —  but  the  method 
is  worthy  of  note  because  it  seems  to  be  the  only  method  currently  employed  which  can 
do  transient  analysis  with  integral  constitutive  equations. 

COMPUTATION  WITH  A  SINGLE-INTEGRAL  MODEL 

It  seems  appropriate  to  present  some  numerical  results  in  this  paper,  but  in  so  doing 
we  must  of  necessity  limit  our  consideration  to  those  results  with  which  the  author  is 
most  familiar.  This  means  that  we  will  look  at  the  extreme  right-wing  method  of  Table  1, 
applied  to  a  flow  which  is  essentially  a  measurement  flow.  We  first  turn  our  attention  to  an 
overview  of  some  of  the  problems  associated  with  the  implementation  of  a  particle-tracking, 
pure  v,  p  method.  More  details  may  be  found  in  refs.  7  and  10. 

Evaluating  the  History- Integral 

The  most  basic  problem  to  be  dealt  with  in  devising  a  numerical  method  for  integral 
constitutive  equations  is  the  numerical  approximation  of  the  history-integrals  in  the  stress 
calculator,  eq.  (4).  Consider  a  typical  memory  function,  m(r);  the  most  efficient  way  to 
compute  the  history-integral,  presuming  the  components  of  the  strain  measure  are  known, 
is  via  a  Gaussian  formula, 

N,. 

S0(r  )m(r )  dr  as  u>*S0 (r*)  (11) 

k-\ 

The  formula  is  in  the  form  of  a  weighted  sum  with  weights,  u>*.  of  the  integrand  evaluated 


The  middle-of-the-road  methods  in  Table  1.  at  first  sight,  seem  applicable  only  to 
differential  models  and  thus  seem  to  be  left-wingers  at  heart,  but  on  closer  examination, 
they  point  the  way  to  a  unification  of  the  two  extremes.  The  middle-of-the-road  methods 
are  currently  under  development  by  Tanner  and  coworkers  [18],  and  Variant  2  currently 
exists  only  as  a  boundary-integral  method,  but  its  obvious  extension  to  a  finite  element 
method  make  it  appropriate  to  this  discussion.  Both  Variants  1  and  2  begin  with  an 
estimated  velocity  field  and  solve  alternate  hyperbolic  and  elliptic  problems  for  stress  and 
velocity,  using  ideas  related  to  the  obsevation  of  eqs.  (10a).  Variant  1  uses  a  discrete 
hyperbolic  method  to  solve  for  the  stress,  and  Variant  2  uses  the  fact  that,  in  the  current 
velocity  iterate,  the  streamlines  are  characteristics,  and  thus  the  stress  equations  can  be 
reduced  to  ordinary  differential  equations  along  them.  These  methods  do  not  alleviate 
the  drawbacks  of  the  extreme  methods,  but  rather  represent  a  different  balance  of  trade¬ 
offs.  They  require  separate  extra  stress  calculation,  but  fewer  unkowns  in  each  phase 
than  the  left-wing  methods.  Each  added  extra  stress  requires  an  extra  pass  through  the 
stress-equation  solver,  not  more  simultaneously  active  unknowns.  The  middle-of-the-road 
methods  require  hyperbolic  data,  but  Newton’s  method  and  transient  analysis  do  not  seem 
as  forbidding  as  with  right-wing  methods. 

However,  the  most  important  aspect  of  the  middle-of-the-road  methods  is  the  bridge 
between  the  extreme  methods  they  illuminate.  Combination  of  Variant  2  and  the  ana¬ 
lytic  streamline  version  of  the  right-wing  method  seem  to  show  real  promise  of  allowing 
the  development  of  a  “super method,”  applicable  to  either  differential  or  integral  models, 
which  shares  as  much  code  as  possible  in  the  computations  with  either  kind  of  constitutive 
equation.  The  common  portion  of  the  code  constructs  the  mesh  and  guides  the  nonlinear 
iterations  (among  other  functions)  and  constructs  streamlines  by  the  procedures  already  in 
place  in  the  analytic  streamline  version  of  the  right-wing  methods.  Depending  on  whether 
the  constitutive  equation  is  integral  or  differential,  a  stress-integral  evaluator  or  a  stress 
ODE  integrator  is  plugged  in  to  evaluate  the  stress  along  the  streamlines.  The  super¬ 
method  will  take  some  time  to  develop,  since  solving  hyperbolic  equations  numerically  is 
not  a  trivial  task.  For  example,  one  can  forsee  some  difficulty  associated  with  the  spec¬ 
ification  of  initial/end  conditions  for  the  stress  ODEs  on  the  streamlines.  For  incoming 
streamlines,  stress  can  probably  can  be  specified  as  data,  but  on  closed  streamlines,  peri¬ 
odic  end-conditions  seem  appropriate.  How  to  do  this  without  excessive  “shooting”  is  not 
clear  at  present.  But  the  possibility  of  the  existence  of  a  supermethod  points  out  that  the 


slot  than  they  are  in  Stokes-flow.  Careful  comparison  of  Figures  4,  6,  and  7  show  that  the 
streamline  nearest  the  slot  is  in  virtually  the  same  position  in  both  Stokes-flow  and  the 
non-Newtonian  flow,  but  that  in  both  Figures  6  and  7,  the  streamline  dips  measureably 
father  down.  Also,  in  distinct  contrast  to  the  flow  of  Figure  5,  numerical  results  reported 
elsewhere  [4,10],  and  flow  visualizations  [4],  the  vortices  in  Figures  6  and  7  are  higher  up 
and  nearer  the  slot  mouth  than  they  are  in  Stokes-flow.  It  seems  that  with  the  concocted 
constitutive  equation,  the  flow  in  the  slot  interacts  with  and  disturbs  the  channel-flow  to 
a  much  greater  extent  than  it  does  in  the  Curtiss-Bird  case. 

It  would  be  easy  to  dismiss  the  differences  between  the  two  constitutive  equations  dis¬ 
cussed  here  as  an  unphysical  consequence  of  the  artificially  concocted  constitutive  equation; 
however,  behavior  like  that  observed  with  the  concocted  equation  seems  to  be  character¬ 
istic  of  all  the  integral  constitutive  equations  the  author  has  tested,  except  the  Curtiss- 
Bird.  And  the  behavior  seems  to  be  worse  for  other  equations;  the  tinkering  applied  to 
the  Johnson-Segalman  model  has  produced  the  second  best  constitutive  equation  to  the 
Cirtiss-Bird,  measured  by  the  ability  to  get  to  high  De  and/or  Lp.  The  difference  between 
the  two  constitutive  equations  seems  to  devolve  from  the  amount  of  shear-thinning  at  a 
given  elasticity;  this  can  be  quantified  by  observing  L0  as  a  function  of  o  in  simple  shear. 
The  Curtiss-Bird  equation  leads  to  an  L0  vs.  o  curve  with  positive  curvature,  while  the 
concocted  equation  has  negative  curvature.  This  leads  to  a  prediction  [22],  which  is  ob¬ 
served  numerically,  that  Pe/Ni  should  be  larger  at  fixed  a  for  the  concocted  equation  than 
it  is  for  Curtiss-Bird.  Referring  to  Figures  5  and  7,  we  can  see  that  the  component  of 
cross-stream  normal  force  acting  along  the  streamline  should  be  less  for  the  Curtiss-Bird 
fluid  than  it  is  for  the  concocted  equation,  since  in  the  Curtiss-Bird  case,  the  streamlines 
are  much  more  nearly  parallel  to  the  wall.  Thus  we  are  brought  back  to  the  old  arguement 
of  Tanner  and  Pipkin  [26],  in  which  the  non-zero  values  of  Pe  observed  experimentally 
are  attributed  to  elastic  restoring  forces  acting  on  a  control  volume  over  the  slot.  This 
kind  of  explanation  of  viscoelastic  hole-pressure  seems  to  be  borne  out  in  these  numerical 
results,  at  least  on  a  qualitative  level.  It  is  also  interesting  to  observe  that  the  same  inter¬ 
action  between  slot  and  channel-flow  which  causes  relatively  higher  hole-pressure  seems  to 
cause  numerical  difficulty.  It  is  tempting  to  speculate  that  the  numerical  difficulty  may  be 
due  to  an  impending  breakdown  of  steady  channel-flow  near  the  slot  mouth,  due  to  this 
interaction. 


The  author  believes  (hat  the  numerical  methods  decribed  here  represent  an  auspi¬ 
cious  beginning  to  the  task  of  numerically  modelling  non-Newtonian  flows;  however,  they 
represent  only  a  beginning.  There  needs  to  be  a  better  theoretical  understanding  of  the 
nature  of  the  equations  and  the  resulting  approximation  schemes  before  definitive  progress 
can  be  made.  Recent  work  illuminating  the  nature  of  the  characteristics  of  the  equations 
moves  in  that  direction  and  serves  to  illustrate  common  ground  between  the  different 
classes  of  constitutive  theories  proposed  by  rheologists.  Jt  is  hoped  that  this  can  lead  to  a 
more  unified  numerical  approach  to  this  challenging  and  demanding  class  of  problems.  In 
the  meantime,  current  numerical  techniques,  though  only  partially  understood,  seem  to  be 
capable  of  uncovering  physically  important  consequences  of  non-Newtonian  fluid  behavior. 

Acknowledgements:  The  numerical  techniques  for  integral  constitutive  equations  described 
here  were  developed  jointly  by  the  author  and  B.  Bernstein  (Dept.  Mathematics.  I.  I.  T.). 
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